﻿package com.moody.util;

public class NormalDistributionUtil {
	
	public static final double MAX = 3.09;
	public static final double MIN = 0.01;
	
	/**
	 * 正太分布反函数
	 * @param x
	 * @param max
	 * @param min
	 * @param mid
	 * @param deepth
	 * @return
	 */
	public static double inverse_Ni(double x){
		double Y = -1 * Math.log(4 * x * (1-x));
		double X1 = Math.sqrt(Y *(2.0611786 - 5.7262204/(Y + 11.640595)));
		double A = 0;
		if(X1 > 4.5) A = -0.000637 * X1*X1 - 0.010437*X1 + 0.059374;
		double X0 = X1 + A;
		double U0 = Math.exp(-1*X0/2)/ Math.sqrt(2*Math.PI);
		return X0 - ((Ni(X0)-x)/U0)*(1-X0*(Ni(X0)-x)/(2*U0));
	}
	public static double NormSDist(double a){
	    double p, b1, b2, b3, b4, b5;
	    
	    p = 0.2316419;
	    b1 = 0.31938153;
	    b2 = -0.356563782;
	    b3 = 1.781477937;
	    b4 = -1.821255978;
	    b5 = 1.330274429;

	    double x = Math.abs(a);
	    double t = 1 / (1 + p * x);
	    
	    double NormSDist = 1 - (1 / (Math.sqrt(2 * Math.PI)) * Math.exp(-a * a / 2)) * (b1 * t + b2 * t * t + b3 * Math.pow(t, 3) + b4 * Math.pow(t, 4) + b5 * Math.pow(t, 5));
	    
	    if (a < 0){
	    	 NormSDist = 1 - NormSDist;
	    }
	    return NormSDist;
	}
	public static double NormSInv(double p){
		double normSInv;
		double a1 = -39.6968302866538, a2 = 220.946098424521, a3 = -275.928510446969;
	    double a4 = 138.357751867269, a5 = -30.6647980661472, a6 = 2.50662827745924;
	    double b1 = -54.4760987982241, b2 = 161.585836858041, b3 = -155.698979859887;
	    double b4 = 66.8013118877197, b5 = -13.2806815528857, c1 = -7.78489400243029E-03;
	    double c2 = -0.322396458041136, c3 = -2.40075827716184, c4 = -2.54973253934373;
	    double c5 = 4.37466414146497, c6 = 2.93816398269878, d1 = 7.78469570904146E-03;
	    double d2 = 0.32246712907004, d3 = 2.445134137143, d4 = 3.75440866190742;
	    double p_low = 0.02425, p_high = 1 - p_low;
	   
	    double q,r;
	   
	    if(p < 0 || p > 1)
	        return 0;
	    if (p < p_low){
	        q = Math.sqrt(-2 * Math.log(p));
	        normSInv = (((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
	    }
	    else if (p <= p_high){
	        q = p - 0.5;
	        r = q * q;
	        normSInv = (((((a1 * r + a2) * r + a3) * r + a4) * r + a5) * r + a6) * q / (((((b1 * r + b2) * r + b3) * r + b4) * r + b5) * r + 1);
	    }
	    else{
	        q = Math.sqrt(-2 * Math.log(1-p)); 
	        normSInv = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
	    }
	    return normSInv;
	}
	

	/**
	 * 正态分布函数近似值，使用6项近似余函数
	 */
	private static double Fi_erf_6(double x) {
		double a = Math.abs(x);
		return 0.5 * (1 + erf_6(a / Math.sqrt(2)));
	}

	/**
	 * 正态分布函数六项级数近似余函数
	 */
	private static double erf_6(double x) {
		double a[] = { 0.070523084, 0.0422820123, 0.0092705272, 0.0001520143,
				0.0002765672, 0.0000430638 };
		double t = 0;
		for (int i = 0; i < 6; i++) {
			t = t + a[i] * Math.pow(x, i + 1);
		}
		return 1 - Math.pow(1 + t, -16);
	}

	/**
	 * 正态分布函数值
	 */
	public static double Ni(double x) {
		return x == 0 ? 0.5 : (x > 0 ? Fi_erf_6(x) : 1 - Fi_erf_6(x));
	}
}
